In [ ]:
# STAT 415/615 Regression (M. Baron)
# Python Lab 1. Introduction and Regression Examples
# Example 1. US population.
In [2]:
# Set a working directory. The name of yours will be different from mine. Then read the CSV (comma-separated values) file.
import os
os.chdir(r"C:\Users\baron\Documents\Teach\615 Regression\Data")
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
In [4]:
P = pd.read_csv("USpop.csv")
P
Out[4]:
| Year | Population | |
|---|---|---|
| 0 | 1790 | 3.9 |
| 1 | 1800 | 5.3 |
| 2 | 1810 | 7.2 |
| 3 | 1820 | 9.6 |
| 4 | 1830 | 12.9 |
| 5 | 1840 | 17.1 |
| 6 | 1850 | 23.2 |
| 7 | 1860 | 31.4 |
| 8 | 1870 | 38.6 |
| 9 | 1880 | 50.2 |
| 10 | 1890 | 63.0 |
| 11 | 1900 | 76.2 |
| 12 | 1910 | 92.2 |
| 13 | 1920 | 106.0 |
| 14 | 1930 | 123.2 |
| 15 | 1940 | 132.2 |
| 16 | 1950 | 151.3 |
| 17 | 1960 | 179.3 |
| 18 | 1970 | 203.3 |
| 19 | 1980 | 226.5 |
| 20 | 1990 | 248.7 |
| 21 | 2000 | 281.4 |
| 22 | 2010 | 308.7 |
| 23 | 2020 | 331.5 |
In [5]:
plt.scatter(P["Year"], P["Population"])
plt.xlabel("Year")
plt.ylabel("Population")
plt.show()
In [7]:
# This is called a scatterplot. Now fit a linear regression line, the closest to the given data points among all straight lines. Then draw this line on our scatterplot and predict the US population for year 2030 using this linear model.
X = sm.add_constant(P["Year"])
y = P["Population"]
regr = sm.OLS(y, X).fit()
plt.scatter(P["Year"], P["Population"])
plt.plot(P["Year"], regr.predict(X), linewidth=3)
plt.show()
In [10]:
new_data = pd.DataFrame({"Year": [2030]})
new_X = pd.DataFrame({
"const": 1,
"Year": new_data["Year"]
})
regr.predict(new_X)
Out[10]:
0 291.534058 dtype: float64
In [11]:
# Clearly, linear regression is not the best choice here because apparently, the population does not grow linearly.
# Our linear model underestimates population in 1800s and 2000s and overestimates in the middle of the graph. As a result,
# the prediction of 291.5 million people two years from now is ridiculous because during the latest US Census in 2020,
# the US population was estimated to be 331.5 million.
# At the same time, we see from the summary of our model (below) that both the slope and the intercept of our regression
# line are significant, and overall, the model explained 92.14% of the total variation, which is generally considered rather good.
print(regr.summary())
OLS Regression Results
==============================================================================
Dep. Variable: Population R-squared: 0.921
Model: OLS Adj. R-squared: 0.918
Method: Least Squares F-statistic: 257.8
Date: Sat, 29 Aug 2026 Prob (F-statistic): 1.24e-13
Time: 16:46:11 Log-Likelihood: -114.71
No. Observations: 24 AIC: 233.4
Df Residuals: 22 BIC: 235.8
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const -2600.4834 169.142 -15.375 0.000 -2951.262 -2249.705
Year 1.4246 0.089 16.056 0.000 1.241 1.609
==============================================================================
Omnibus: 3.914 Durbin-Watson: 0.096
Prob(Omnibus): 0.141 Jarque-Bera (JB): 2.503
Skew: 0.597 Prob(JB): 0.286
Kurtosis: 1.963 Cond. No. 5.25e+04
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 5.25e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
In [13]:
# Let’s try a more advanced quadratic model. Among all quadratic models (that is, all parabolas in the world),
# this one is the closest to our data points. Then, we calculate predicted values “Yhat” that are obtained by
# plugging all years 1790, 1800, …, 2020 into the resulting quadratic polynomial and plot this curve in blue.
X_quad = pd.DataFrame({
"Year": P["Year"],
"Year2": P["Year"] ** 2
})
X_quad = sm.add_constant(X_quad)
quad = sm.OLS(P["Population"], X_quad).fit()
print(quad.summary())
OLS Regression Results
==============================================================================
Dep. Variable: Population R-squared: 0.999
Model: OLS Adj. R-squared: 0.999
Method: Least Squares F-statistic: 1.392e+04
Date: Sat, 29 Aug 2026 Prob (F-statistic): 1.62e-33
Time: 16:48:36 Log-Likelihood: -58.939
No. Observations: 24 AIC: 123.9
Df Residuals: 21 BIC: 127.4
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 2.171e+04 522.079 41.579 0.000 2.06e+04 2.28e+04
Year -24.1296 0.549 -43.982 0.000 -25.270 -22.989
Year2 0.0067 0.000 46.585 0.000 0.006 0.007
==============================================================================
Omnibus: 10.607 Durbin-Watson: 1.280
Prob(Omnibus): 0.005 Jarque-Bera (JB): 8.697
Skew: -1.227 Prob(JB): 0.0129
Kurtosis: 4.635 Cond. No. 3.09e+09
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 3.09e+09. This might indicate that there are
strong multicollinearity or other numerical problems.
In [14]:
Yhat = quad.predict(X_quad)
plt.scatter(P["Year"], P["Population"])
plt.plot(P["Year"], Yhat, linewidth=3)
plt.show()
In [15]:
# Wow, it fits the data really well! Also, now we explained 99.9% of the total variation. Only twice, during
# the WWII and right after, the US population appeared just a little below the curve. However, the “baby boom”
# came, and it quickly caught up with the quadratic trend.
# Now, use the new model to make forecasts:
new_data = pd.DataFrame({"Year": [2030, 2040, 2050]})
new_X = pd.DataFrame({
"Year": new_data["Year"],
"Year2": new_data["Year"] ** 2
})
new_X = sm.add_constant(new_X)
quad.predict(new_X)
Out[15]:
0 364.194713 1 395.879662 2 428.906038 dtype: float64
In [ ]:
# Our quadratic model predicts the US population of 364 mln. people in year 2030, 396 mln. in 2040, and 429 mln. in 2050.
# Let’s wait and see. The current official estimate for 2026 is 349 mln.
In [16]:
# Example 2. US Presidents.
# The given data set contains Presidents from Andrew Johnson to Lyndon Johnson. It compares the actual number of years
# each President lived after his first inauguration against the average number of years that a man of the same age would
# live during the same time.
Pres = pd.read_csv("presidents.csv")
Pres
Out[16]:
| name | expected | actual | |
|---|---|---|---|
| 0 | ANDREW JOHNSON | 17.2 | 10.3 |
| 1 | ULYSSES S. GRANT | 22.8 | 16.4 |
| 2 | RUTHERFORD B. HAYES | 18.0 | 15.9 |
| 3 | JAMES A. GARFIELD | 21.2 | 0.5 |
| 4 | CHESTER A. ARTHUR | 20.1 | 5.2 |
| 5 | GROVER CLEVELAND | 22.1 | 23.3 |
| 6 | BENJAMIN HARRISON | 17.2 | 12.0 |
| 7 | WILLIAM MCKINLEY | 18.2 | 4.5 |
| 8 | THEODORE ROOSEVELT | 26.1 | 17.3 |
| 9 | WILLIAM H. TAFT | 20.3 | 21.2 |
| 10 | WOODROW WILSON | 17.1 | 10.9 |
| 11 | WARREN G. HARDING | 18.1 | 2.4 |
| 12 | CALVIN COOLIDGE | 21.4 | 9.4 |
| 13 | HERBERT C. HOOVER | 19.0 | 35.6 |
| 14 | FRANKLIN D. ROOSEVELT | 21.7 | 12.1 |
| 15 | HARRY S. TRUMAN | 15.3 | 27.7 |
| 16 | DWIGHT D. EISENHOWER | 14.7 | 16.2 |
| 17 | JOHN F. KENNEDY | 28.5 | 2.8 |
| 18 | LYNDON B. JOHNSON | 19.3 | 9.2 |
In [17]:
plt.scatter(Pres["expected"], Pres["actual"])
plt.xlabel("expected")
plt.ylabel("actual")
plt.show()
In [18]:
# If presidents are just ordinary people, and the fact of presidency does not affect their lifetime,
# then the points should be close to the bisector line y = x. But surprisingly, here we see that this
# is not the case. OK, let’s fit a linear regression line.
X = sm.add_constant(Pres["expected"])
y = Pres["actual"]
reg = sm.OLS(y, X).fit()
plt.scatter(Pres["expected"], Pres["actual"])
plt.plot(Pres["expected"], reg.predict(X), linewidth=3)
plt.show()
In [19]:
print(reg.params)
const 23.386720 expected -0.506074 dtype: float64
In [ ]:
# What do we see? The slope is negative! It means that presidents who were expected
# to live longer actually passed away sooner. Perhaps, it is because three of the
# listed presidents were assassinated – Garfield, McKinley, and Kennedy (#4, 8, 18 in the data set).
# One may argue that they should be excluded from modeling (although others may disagree with this
# referring to the US history - the risk of being assassinated is much higher when you become
# a US President, so these three may be a part of the trend and not an outlier). But the whole
# picture is not that different even without the assassinated presidents.
In [20]:
# R uses 1-based row numbers; Python uses 0-based row numbers. For example, President #4 is in row #3 in Python.
Z = [3, 7, 17]
Pres_no_assassinations = Pres.drop(index=Z)
X = sm.add_constant(Pres_no_assassinations["expected"])
y = Pres_no_assassinations["actual"]
reg = sm.OLS(y, X).fit()
print(reg.summary())
# The slope is still negative. We can see from a high p-value actually, it is not significant, based on our data.
# That is, there is no evidence that the slope is not 0. In other words, knowing the age of a US President at his
# or her first inauguration and the life expectancy of a person of the same age does not really help to predict
# the life expectancy of the President. A surprising conclusion, I think.
OLS Regression Results
==============================================================================
Dep. Variable: actual R-squared: 0.000
Model: OLS Adj. R-squared: -0.071
Method: Least Squares F-statistic: 1.994e-05
Date: Sat, 29 Aug 2026 Prob (F-statistic): 0.997
Time: 16:54:50 Log-Likelihood: -56.378
No. Observations: 16 AIC: 116.8
Df Residuals: 14 BIC: 118.3
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 15.3849 14.970 1.028 0.322 -16.723 47.493
expected -0.0034 0.763 -0.004 0.997 -1.641 1.634
==============================================================================
Omnibus: 3.496 Durbin-Watson: 2.359
Prob(Omnibus): 0.174 Jarque-Bera (JB): 1.759
Skew: 0.797 Prob(JB): 0.415
Kurtosis: 3.312 Cond. No. 134.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
C:\Users\baron\AppData\Local\anaconda3\Lib\site-packages\scipy\stats\_axis_nan_policy.py:531: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16 res = hypotest_fun_out(*samples, **kwds)
In [ ]:
# (Note that the pink comment is not an error. It is just a warning of a sample too small for an accurate kurtosis testing)